Vacancy diffusion in colloidal crystals as determined by dynamical density functional 

theory and the phase field crystal model 
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A two-dimensional crystal of repulsive dipolar particles is studied in the vicinity of its melting 
transition by using Brownian dynamics computer simulation, dynamical density functional theory 
and phase-field crystal modelling. A vacancy is created by taking out a particle from an equilibrated 
crystal and the relaxation dynamics of the vacancy is followed by monitoring the time-dependent 
one-particle density. We find that the vacancy is quickly filled up by diffusive hopping of neighbour- 
ing particles towards the vacancy center. We examine the temperature dependence of the diffusion 
constant and find that it decreases with decreasing temperature in the simulations. This trend is 
reproduced by the dynamical density functional theory. Conversely, the phase field crystal calcula- 
tions predict the opposite trend. Therefore, the phase-field model needs a temperature-dependent 
expression for the mobility to predict trends correctly. 
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I. INTRODUCTION 

Most of the mechanical properties of crystals depend 
crucially on the presence of crystalline defects. For ma- 
terial processing it is therefore of principal importance 
to understand and control the defect concentration and 
dynamics. The nature and dynamics of defects is much 
easier to classify for crystalline sheets in two spatial di- 
mensions. In this case, it is known for a long time, that 
the formation and unbinding of topological defects pro- 
vides an efficient way of melting according to the two- 
stage scenario proposed by Kosterlitz-Thouless-Nelson- 
Halperin- Young (KTNHY) [l|. Defects can also be accu- 
mulated near edges of crystalline sheets and do occur for 
two-dimensional crystals on more complicated manifolds 
ii. 

Defects are highly dynamic: Whereas the structure of 
a crystal is static over long time scales, defects undergo 
diffusion in the crystalline background. The diffusive 
dynamics of individual point defects were observed di- 
rectly in two-dimensional colloidal suspensions of charged 
micro-spheres by video microscopy [j, [sj and they were 
confirmed and further analyzed by computer simula- 
tions [Hi- The dynamics of defects were also explored 
by real-space methods in a two-dimensional crystal of 
weakly damped dust particles in a plasma by Nosenko 
and coworkers 041], see also [loj . 

Describing the defect dynamics by a microscopic the- 
ory is still a formidable challenge, in particular close to 
melting. Such a theory should contain the solid and fluid 
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phase and give a reliable picture on the defect concen- 
tration and its dynamics. Recent progress in this respect 
has been made by classical density functional approaches 
of freezing |ll| - [l4| . For Brownian particles, the density 
functional approach can be generalized towards dynam- 
ics [Tsl - flTj and the dynamics of solidification has been 
examined in two dimensions [Tsl . [Tgj . Similar in line, a 
more coarse-grained phase-field-crystal model has been 
proposed to describe crystal growth [20l - l23j and the de- 
fect structure and dynamics, for various applications see 
e.g. [24l - ]33 |. However a systematic exploration of defect 
dynamics by such a density functional theory has not yet 
been performed nor has the reliability of the phase-field- 
crystal model been systematically checked as far as the 
trends of defect dynamics are concerned. 

Here, we study the dynamics of vacancies in a two- 
dimensional colloidal crystal by using Brownian dynam- 
ics computer simulations, dynamical density functional 
theory and the phase-field crystal approach and thereby 
test the ability of the theoretical approaches to quali- 
tatively reproduce the observations made in the simula- 
tions. The model system we use here is a two-dimensional 
suspension of dipolar colloids. This system has been re- 
alized experimentally as superparamagnetic particles at 
an air- water interface (35| . When exposing superparam- 
agnetic particles to an external magnetic field perpendic- 
ular to the plane, their induced magnetic dipole moment 
leads to an effective repulsive interaction whose ampli- 
tude can be tuned by the magnetic field strength. At suf- 
ficiently high field strength, the system crystallizes into a 
two-dimensional triangular (i.e. hexagonal) crystal. This 
system has been studied extensively by computer sim- 
ulations and by the aforementioned dynamical de nsity 
functional theory [l^ and phase field crystal models fl9| . 

Out of a perfect triangular crystal, some particles are 
removed and the relaxation of the resulting defect and 
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its mobility are extracted by monitoring the one-particle 
density function of time. We confirm by simula- 
tion that the defect mobility is increasing with increasing 
temperature, as was already observed for charged par- 
ticles by Li'bal et al Q. This behavior is qualitatively 
and semi-quantitatively reproduced in dynamical den- 
sity functional theory based on a static Ramakrishnan- 
Yussouff-like density functional [36i, |37j] . The phase-field- 
crystal model, on the other hand, fails to predict the 
trend of the temperature-dependence of the mobilities. 
This is mainly attributed to the constant kinetic pref- 
actor involved in the phase-field-crystal approach. For 
predicting the trend, a temperature-dependent corrective 
mobility is needed for the phase-field-crystal model. 

This paper is organized as follows: in section II, we 
briefly propose the different approaches used in this pa- 
per. Results are presented in section III and we conclude 
in section IV. 



II. THEORETICAL MODELS 

Dynamical density functional theory and the more 
coarse-grained phase-field-crystal model describe the 
overdamped Brownian dynamics in terms of a continuity 
equation for the deterministic, time-dependent, and en- 
semble averaged one-particle density p(r, t). Note that in 
many applications of the phase-field-crystal model p(r, t) 
is interpreted as a fluctuating density fleld that changes 
for different realizations of the dynamical evolution even 
under the same initial conditions. Here, we will not take 
this approach but regard p(r, t) as a purely deterministic 
quantity. For a more thorough discussion of the physical 
interpretation of p(r, t) in the phase-field-crystal model 
see [ii,!!^. 

The temporal evolution of the density field according 
to dynamical density functional theory [l^ [s^ is given 
by 
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with D the single-particle diffusion constant and fc^T the 
thermal energy. The Helmholtz free energy functional 
Fjjo(r)] is provided by classical density functional theory 
In the crystal, the driving current pW{SF/Sp) obvi- 
ously assumes the hexagonal symmetry of the underlying 
crystal. In the interstitial regions, this current is small 
since the density itself is small. 

Note that Eq. ([T]) can be derived from first principles 
[l^ll^, i.e., from the microscopic Langevin equations of 
motion or from the Smoluchowski equation for the time 
evolution of their respective probability distribution (for 
a review see [l9|). Here, we employ the same approxima- 
tion to the density functional theory of Ramakrishnan 
and Yussouff [l^ already introduced in Refs. [l^ fisj . 



Eq. ([T]) then reads 



p{r, t)^D\ V2p(r, t) + {kBT)-^W ■ [p{r, t)WV{r, t)] 



p(r,i)V / drV(r')4'^(|r-r'|;p) 
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where V{r,t) is the time-dependent external potential. 

(2) 

Fcx{p) and Cq (r; p) arc the excess free energy and the 
direct correlation function of the reference fluid of density 
p, respectively. 

In this work we consider a two-dimensional system of 
magnetic dipoles that are oriented perpendicular to the 
2D plane. The pair potential of two dipoles in the plane 
is given by 



u{r) = uo/r^ 



(3) 



where uq is the interaction strength. The thermodynam- 
ics and structure depend only on one dimensionless cou- 
pling parameter F = u^p^/"^ /ksT ^ where p is the average 
one-particle density and /csT is the thermal energy. 

The two-particle direct correlation function of the fluid 
c^q\y) [l^l has been obtained for a large range of cou- 
pling constants < F < 62.5 from liquid-state integral 
equation theory as previously described [l^ [13, Ell ■ 

In order to measure the diffusion of defects, Eq. ((H is 
numerically solved on a rectangular periodic box of a fine 
grid with ^ 64 grid points per nearest-neighbor distance 
a. A finite difference method with variable time step is 
applied. The convolution integrals are solved using the 
method of fast Fourier transform. 

For the more coarse-grained phase-ficld-crystal model 
we employ the two different versions termed PFCl and 
PFC2 in Ref. [H. The PFCl model constitutes an ap- 
proximation of dynamical density functional theory, as 
introduced above. The last term in Eq. ^ is replaced 
by its gradient expansion. The dynamical equation then 
reads 



p(r, t) = DV2p(r, t) + J p(r, i)V (fesT)" V(r, t) 



(4) 



The parameters Co, C2, and C4 are the fit parameters of 
a parabola to the first maximum of the Fourier transform 
of the two-particle correlation function Cg (fc). The co- 
efficient a = 1.15 is artificially introduced to match the 
melting points of PEG and DDFT. 

The second, more coarse-grained model termed PFC2 
in Ref. [l^, which is frequently used in the phase-field- 
crystal literature, can be obtained from dynamical den- 
sity functional theory by assuming a constant mobility. 
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FIG. 1: Sketch of the setups in the simulation (left) and the- 
ory (right): In the computer simulation we remove one parti- 
cle out of a perfect hexagonal lattice and follow the position 
of the vacancy over time. In dynamical density functional 
theory and the phase-field crystal model we study a quasi- 
one- dimensional relaxation of a depleted central density peaks 
(grey) being replenished by influx of probability density from 
surrounding density peaks. 



/o(r, t) = p in front of the functional derivative in Eq. ([T]) 
and a gradient expansion. The model equation then reads 



+ {kBT)-^V{r, t) - ap (do - CaV' + QV^ ) <^(r, t) 



FIG. 2: a) A vacancy typically appears as one of four dif- 
ferent configurations of particles with 5, 7 or 8 neighbors, 
indicated in light, dark grey and black, respectively, b) A 
typical trajectory of a vacancy (similar to results presented in 
Ref. i). 



(5) 




with (p{r,t) = [p{r,t) — p]/p the dimensionless density 
modulation. 

The equation of motion is solved using the finite differ- 
ence method with a semi-implicit time integration [4l| . 



III. SETUP AND RESULTS 



FIG. 3: The mean square displacement for F = 16.6 (black 
dashed line), F = 21.12 (black dotted line), F = 24 (black 
dash-dotted line), and F = 28.8 (black continuous line) ob- 
tained from Brownian Dynamics computer simulations. The 
grey continuous lines are linear fits. Inset: the vacancy dif- 
fusion constant Dv, calculated from the slope of the mean 
square displacement, as a function of F. 



In the following subsections we qualitatively compare 
the temperature dependence of defect diffusion as ob- 
tained by computational Brownian dynamics simulations 
and as predicted by the dynamical density functional the- 
ory and the phase- field-crystal model 2 (PFC2). 

A. Brownian dynamics computer simulation 

As a reference for the theoretical models we use Brow- 
nian dynamics computer simulations [4^ to quantify the 
diffusion of vacancies for different coupling strengths F 
(well above the melting point at Fm ~ 12[43)- Follow- 
ing Libal et al. Q, we equilibrate a perfectly hexagonal 
crystal of = 2500 particles in a rectangular, almost 
square box employing periodic boundary conditions while 
keeping one particle tagged at the origin (see Fig. [1]). 
Subsequently, the particle at the origin is removed and 



the diffusion of the defect is followed over time. In this 
setup the vacancy concentration is typically relatively 
small (of the order 10~^). The position of the defect 
is determined by a Voronoi construction: The vacancy 
appears as one of several configurations of multiple par- 
ticles with more or less than six neighbors (see Fig. 
The center of mass of these particles is considered as the 
position of the vacancy. 

As was already observed for Yukawa particles in 2D by 
Libal et al. Q , we also find that the defect undergoes dif- 
fusion and that the diffusion constant increases with in- 
creasing temperature, corresponding to a decreasing cou- 
pling strength F (Fig. [3]). The diffusion constant of the 
vacancy Dy ranges between 18.9(prB)^^, for F = 16.6, 
and 9{ptb)-\ for F = 28.8. 
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B. Theory 

In dynamical density functional theory and in the 
phase-ficld-crystal models, crystals appear as strongly 
modulated density fields that have the symmetry of the 
corresponding crystal [l^, [l^ . These density fields are 
mechanically and thermodynamically stable at low tem- 
perature or high coupling strength f37| . In an equilib- 
rium density field, the integrated density field over one 
Wigner-Seitz cell is equal to the probability to find a par- 
ticle at the corresponding lattice site. In the approxima- 
tion to the density functional theory by Ramakrishnan 
and Yussouff for magnetic dipoles in 2D described above 
this number is very close to 1. A number smaller than 
1 can be interpreted as a finite probability to find a va- 
cancy, a number larger than one as a finite probability to 
find an interstitial, respectively. 

Whereas we have addressed the short-time relaxation 
dynamics of crystals in a previous paper [Tsj . we arc 
here concerned with the long-time dynamics of vacan- 
cies. For ease of computation and to assess larger time 
scales we thus start with a slightly different setup, as 
compared to the setup in the computer simulations: In- 
stead of introducing a vacancy with probability 1 at the 
center of a large two-dimensional crystal, which consti- 
tutes a problem of cylindrical symmetry, we study the 
relaxation dynamics of the quasi-one-dimensional setup 
sketched in Fig. [TJ Our reference state is an equilib- 
rium crystalline density field in a periodic rectangular 
box of dimensions Lx x Ly = 2a x 64(-\/3/2)a, where a 
is the equilibrium lattice constant. At time t = we 
reduce the integrated density of all density peaks lying 
on an infinite row of neighboring crystal sites along the 
X axis by a small amount of only 3 percent, thus ren- 
dering the problem quasi-onc-dimcnsional (sec Fig. [TJ. 
Specifically, low-amplitude Gaussians with the same cen- 
ter positions and similar width as the equilibrated den- 
sity peaks were subtracted from the density field. Thus, 
the integrated density of each altered peak is smaller by 
a few percent than its equilibrium value. The temporal 
behavior of this setup is expected to be qualitatively sim- 
ilar to that of a cylindrical setup at late times, i.e., once 
the vacancy has diffused a large distance from its initial 
position [561 . Strictly speaking the long-time diffusion 
constant D should be calculated from the self-part of a 
van-Hove function in the presence of interactions (e.g. by 
using the test particle approach to the DDFT |4a]). For 
a low vacancy density, however, the collective and self- 
dynamics are expected to be similar such that we have 
avoided the significant additional effort to implement the 
tcst-particlc approach. 

The outcome of the dynamical density functional the- 
ory is summarized in Fig. 2] and based on the x- averaged 
density field 

Px{y,t) = L-' J dxpir,t). (6) 
Fig. 2^) displays the equilibrium averaged density field 
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FIG. 4: Results of the dynamical density functional theory: 
a) the initial, equilibrium, x-averaged density profiles p"(j/,t) 
for two difi'erent values of F, F = 40 and F — 62.5, the former 
being close to the melting point at F = 36, b) the difference of 
the density profile at time t and the profile at time measured 
at the positions of the peaks in (a) , such that a set of discrete 
data set is obtained which is in monotonic in y, c) the variance 
o-^(f) of a Gauss function fitted to (b) as a function of time, 
d) the effective diffusion constant calculated as D = ay/{2t). 



Px{y) — Px{y,t — 0), i.e., before the introduction of 
vacancies, for two different coupling strength, F = 40 
and F = 62.5, that are close and far from the freezing 
transition at T f fa 36, respectively [53 • The higher 
coupling strength corresponds to higher and more pro- 
nounced density peaks. Removing a fraction of 0.03 par- 
ticles from the row of peaks at y = leads to a restoring 
density current from neighboring particle rows towards 
the origin. This is represented by the difference between 
the x-averaged density fields of the perturbed and the 
unperturbed systems at their original y-positions (Fig. 

Apx^ pliy)~ Px{y,t). (7) 

For small initial perturbations and at long times the en- 
velope of Apx approaches a Gaussian function, which 
broadens over time. The variance (7v(^) plotted in Fig. 
[4];. As expected, cry{t) shows a linear dependence of t 
at long times. The long-time slope translates into a dif- 
fusion constant given by = (Ty/(2t) (Fig. HJl). In 
agreement with the Brownian dynamics simulations the 
diffusion constant is higher for low coupling constant of 
F = 40 (Dv ~ 7(ptb)^^) than for the high coupling con- 
stant of F = 62.5 {{D^ « 4(prB)"^). 

The same setup is studied in the PFCl and PFC2 mod- 
els. The temporal evolution cry(t) of the width of the 
Gaussian envelope function describing the relaxing den- 
sity field is shown for the PFC2 model in Figure [SJ After 
a fast transient relaxation process its dependence is lin- 
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FIG. 5: Results of the PFC2 model: The mean square dis- 
placement CTv for r = 40 (black dashed line), F = 55 (black 
dash-dotted line), and F = 62.5 (black continuous line). The 
grey continuous lines are linear fits to the curves. Inset: the 
vacancy diffusion constant Dv, calculated from the slope of 
the mean square displacement as a function of F. 



ear in time. Remarkably, the corresponding slope -Dy 
which is presented in the inset of Figure [5] is increasing 
for increasing F conversely to what has been found before 
in simulation and DDFT. 

The diffusion constant Dy increases linearly with F 
from 27(/9tb)~^ corresponding to F = 40 to 33(/9tb)^^ 
and 35.7(ptb)~i for F equal to 55 and 62.5. The PFCl 
model gives the same incorrect trend as the PFC2 mod- 
els. Again, after a transient process the system reaches a 
state where the relaxation towards equilibrium is getting 
diffusive but the slope increases with increasing F. The 
physical reason for the incorrect trend in both variants of 
PFC models is that the PFC has a smoothened density 
profile and therefore allows for a quick diffusive current 
of particles from one lattice site to another. In DDFT, 
on the other hand, the full density profile is kept and the 
density decays to very small values in the interstitial re- 
gion between the density peaks. Thus, particle currents 
between lattice sites are strongly reduced. For increas- 
ing coupling F, the interstitial density drops even fur- 
ther down and therefore reduces vacancy diffusion more. 
This effect is not contained in both PFC models. Here 
rather the only remaining trend is set by the increasing 
interactions which leads to larger inter-particle forces and 
therefore accelerate the dynamics of vacancy diffusion. 
Therefore, to account for the correct defect dynamics, 
the temperature dependence of the mobility prcfactor in 
the PFC models needs a proper fitting to reference data. 
Though this reduces the predictive power of the PFC 
model, it may still be useful for a fast explorative numer- 



ical study for various dynamical effects in solids provided 
this fitting is before a priori. 

IV. CONCLUSIONS 

In conclusion, we have used dynamical density func- 
tional theory, phase-field-crystal theory and particle- 
resolved Brownian dynamics computer simulations to 
calculate the diffusion of defects in a two-dimensional 
crystal of repulsive dipoles. The typical diffusion coef- 
ficient of defect motion is expected to decrease with in- 
creasing system temperature as confirmed by the simu- 
lation data. This trend is reproduced in the dynamical 
density functional theory but not in the phase field crys- 
tal calculations. These findings show that the PFC model 
requires a fitting of the kinetic mobilities as a function of 
the thermodynamic parameters if a realistic description 
of trends is required to predict material properties. Given 
the fact that the efficiency of the PFC model is in general 
achieved by an optimal fitting procedure, also for struc- 
tural predictions [l^, such a phenomenological input of 
the mobility is an acceptable fact. However, it shows that 
clearly the dynamical density functional theory is more 
appropriate to predict the microscopic time evolution as 
a first-principle theory for Brownian systems. 

Future work should address the dynamics of dipolar 
mixtures of colloidal particles with different dipole mo- 
ments [4?! where an equilibrium density functional for a 
binary mixture (isj is needed. These mixtures show more 
complex possibilities of mixed crystals as a function of the 
asymmetry in their dipole moments [49| . In this case one 
will expect different diffusion coefficients for different de- 
fects topologies. Moreover, one should do a similar calcu- 
lation for hard discs, for which a very accurate functional 
based on fundamental measure theory (sij was proposed 
recently [soj . A similar comparison can be performed in 
three dimensions, e.g. for hard spheres, where the phase 
field crystal model has been tested against density func- 
tional theory recently [l^ . Finally it would be interesting 
to consider particles with orientational degrees of free- 
dom which form liquid crystals with interesting defect 
structures. In fact, phase-field-crystal models for liquid 
crystals have been developed H, HI] and applied [53 to 
freezing recently which opens the way to describe defect 
dynamics in liquid crystalline states. 
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